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I. INTRODUCTION 

A great deal of the dynamics of open systems can be 
described, to a reasonable degree of accuracy, by Marko- 
vian quantum master equations. Important examples 
are given by the weak-coupling interaction of radiation 
with matter in atomic physics and quantum optics GIB- 
However, non-Markovian quantum dynamics 0, Gl- S 13 
is known to play a significant role in many applications 
of the theory of open quantum systems currently un- 
der discussion in the literature, e.g. the dynamics of the 
atom laser ||, environment-induced decoherence at low 
temperatures (for an example, see 111), and quantum de- 
vices interacting with a spin bath |lCj . 

Quantum Monte Carlo techniques have been shown to 
provide efficient numerical tools for the treatment of the 
dynamics of open systems in the Markovian regime 
In these techniques one constructs a stochastic dynam- 
ics for the open system's state vector ip(t) such that the 
reduced density matrix ps(t) of the open system is re- 
covered through the expression ps(t) = E(|-0(*)) 
where E denotes the expectation value or ensemble aver- 
age of the underlying process. This is the standard Monte 
Carlo wave function method which has been widely used 
in many physical problems of quantum optics and con- 
densed matter theory. 

The idea of the Monte Carlo wave function method can 
be extended to the treatment of non-Markovian quan- 
tum processes which cannot be described by a Markovian 
quantum master equation. One such method 0] is based 
on a stochastic integro-differential equation for the wave 
function involving a non-local retarded memory kernel. 
The solution of non-local equations of motion can be cir- 
cumvented by employing a pair ipi(t), ip2 (t) of random 
wave functions of the open system and by expressing the 
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reduced density matrix with the help of the mean value 
p s (t) = E(|^i(£))(V> 2 (£)|) 0. This method of propagat- 
ing a pair of wave functions requires the construction of 
an appropriate time-local non-Markovian master equa- 
tion. Such an equation can be obtained with the help of 
the time-convolutionless (TCL) projection operator tech- 
nique which leads to a systematic perturbation expansion 
for the time-dependent generator of the master equation. 
However, for strong system-environment couplings calcu- 
lations based on the TCL expansion become extremely 
complicated and the derivation of an appropriate TCL 
generator of high order is, in general, not feasible in prac- 
tice. A further possibility is to use an explicit expression 
for the influence functional of the open system to ob- 
tain stochastic differential equations for a pair of random 
wave functions 0]. This method is, however, restricted 
to Gaussian reservoirs and linear dissipation. 

In this paper the details of a new method proposed 
in 0] are presented, which allows to attack the prob- 
lem of non-Markovian quantum evolution by means of 
a Monte Carlo wave function technique. The basic idea 
is to introduce a pair |$i(t)), |$2(i)) of random states 
of the total system, with the aim of a stochastic formu- 
lation of the exact von Neumann dynamics of the com- 
posite system. A similar idea has been used recently to 
construct an exact diffusion process for a pair of one- 
particle wave functions describing systems of identical 
Bosons 0] and Fermions 0]. Here, the state vector dy- 
namics is assumed to represent a piecewise deterministic 
process (PDP). This is a Markovian jump process with 
smooth, deterministic evolution periods between succes- 
sive jumps. The stochastic states of the total system 
are supposed to be tensor product states of the form 
l^i) = "01 ® Xi an d |<1>2) = ip2 (8 Xi- The method thus 
involves four stochastic state vectors, namely a pair 
-02 of state vectors of the open system, and a pair xi , Xi 
of state vectors of the environment. The open system's 
reduced density matrix can then be represented in terms 
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of the expectation value 

PS (t) = E(\Mt)){Mt)\{x2(t)\xi(t)))' (i) 

Contrary to the standard methods mentioned above, this 
representation employs an average over the product of 
two quantities: The dyadic \ipi) (ip2 | of a pair of state vec- 
tors of the open system, and the scalar product (X2IX1) 
of a corresponding pair of environment states. It will 
be shown that this representation allows to design a 
Markovian stochastic process which unravels the full non- 
Markovian behavior of the reduced density matrix. 

The paper is structured as follows. Section II con- 
tains the general construction of the PDP representing 
the exact von Neumann dynamics of the composite sys- 
tem, an investigation of the dynamics of the fluctuations 
of the stochastic process, as well as a detailed descrip- 
tion of the Monte Carlo algorithm of the open system 
dynamics. The example of the non-perturbative decay of 
a two-state system into a bosonic reservoir is discussed 
in Sec. III. This section contains numerical simulations 
of the non-Markovian dynamics of the decay into a reser- 
voir in the regime of strong couplings and corresponding 
long memory times. The quantum dynamics of a specific 
spin bath model is investigated in Sec. IV. This model 
describes the interaction of a single electron spin in a 
quantum dot with an external magnetic field and a bath 
of nuclear spins. Section V contains the conclusions and 
indicates various potential generalizations of the stochas- 
tic method. 



II. GENERAL FORMULATION OF THE 
METHOD 

A. Construction of the PDP 

We investigate the general situation of an open system 
with underlying Hilbert space Hs, which is coupled to 
an environment with Hilbert space He- The state space 
of the composite, total quantum system is given by the 
tensor product Hs <8> He- Working in the interaction 
picture we write the Hamiltonian describing the system- 
environment interaction as 

H I {t)=Y J A a {t)®B a {t). (2) 

a 

The A a (t) and the B a (t) are interaction picture opera- 
tors acting in Hs and He, respectively. The evolution of 
the density matrix p(t) of the total system is then gov- 
erned by the von Neumann equation (h = 1), 

±p{t) = - l [H I {t),p{t)]. (3) 

Our central goal is to construct a representation of p(t) 
in terms of the expectation value 



which is determined through a pair |$i(i)), | *&2 ) ) of 
stochastic state vectors of the composite quantum sys- 
tem. Equivalently, one may define the quantity 

R(t) = \*i(t)){* 3 (t)\, (5) 

which is a random operator on Hs ®He, and write the 
density matrix as the mean value of this operator, that 
is p(t) =E(R(t)). 

In the following we suppose that the stochastic state 
vectors |3>„(i)) (y = 1,2) introduced in Eq. are di- 
rect products of certain system states ipv(i) G Hs and 
environment states Xv(fy € He, that is we have 

|$„(t)> =Mt)®X»(t), v=l,2. (6) 

The reduced density matrix ps (t) of the open system is 
defined through the partial trace over the variables of the 
environment, ps(t) — tr^p(t). In view of Eqs. (@J and 
© this definition immediately leads to the relation 

It is important to realize that a representation of the 
form given in Eqs. I@J and © is possible for any ini- 
tial state p(t = 0). This means that any given density 
matrix p = p(0) of the composite quantum system can 
be written as the mean value p = E(|$i)($2|), hi which 
the random states \$> v ) are direct products of the form 
In particular, it is not necessary to demand that 
p describes an initial state without system-environment 
correlations. 

A formal proof of this statement may be carried out 
as follows. One first observes that a sequence of pairs 
(|3>i), | $2)) °f state vectors, which occur with corre- 
sponding probabilities p\, gives rise to the expectation 
value 

/0 = E(|$ 1 )($ 2 |)=^ W |^)(^|. (7) 

A 

Of course, p\ provides a probability distribution satis- 
fying p\ > and J2\P^ = 1- Introducing new states 
through the relation \^f^) = we can write 

p = X;i^X*al- (8) 

A 

Thus, to prove the above statement we have to show that 
any given density matrix p of the composite quantum 
system can be brought into the form (JSJ), whereby the 
I**) must be direct products. To demonstrate that this 
is in fact possible we introduce an ON-basis {ipi} in Hs 
and an ON-basis {xn} in He and write the given p as 
follows, 

P = X! Pijnm\lpi) {i>j \®\Xn) (Xrn \ , (9) 
ijnm 

where 



p(t) = E(|* 1 (i)>(*a(*)|), (4) 



Pijnm = {4>iXn\p\4>jXm) = \Pijnm\e 2% ' 
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Next, one introduces a collective index A = (ijnm) and 
defines the states 

1*1 ) = \/\pii nm \e +iVi '™il>i®Xn, (10) 

l*2> = \J\pijnm\e- i!Piinm ll>j®Xm, (H) 

which allow one to write Eq. JJjJ m t ne desired form JSJl. 
This completes the proof since the states (|10|) and <|11[) 
are indeed direct products. 

The aim is now to construct an appropriate stochas- 
tic process for the state vectors | $„(£)) which exactly 
reproduces the von Neumann equation @ through the 
expectation value Q). As mentioned in the Introduction 
we suppose that the time-evolution represents a piecewise 
deterministic process (PDP). A convenient way of formu- 
lating a PDP is to write stochastic differential equations 
for the random variables. The foundations of the calcu- 
lus of PDPs and its applications to the quantum theory 
of open systems may be found in • In view of the rep- 
resentation © the stochastic dynamics can be defined 
in terms of stochastic differential equations for the state 
vectors and Xv(t), 

dip v (t) = F u dt + dJ u , (12) 
d Xv (t) = G„dt + dK u . (13) 

These equations reflect the general structure of a PDP: 
The terms F v dt and G v dt represent the deterministic evo- 
lution periods, the drift of the process, while the terms 
dJ v and dK v provide the contributions from the random, 
instantaneous jumps of the process. These jump contri- 
butions are taken to be of the form 

dJ v = J2(- iL ^A a -I)i/; v dN au (t), (14) 

a 

dK v = J2 {M a „B a - I) X vdN ay {t). (15) 

a 

Here, / denotes the identity operator and L a „, M av 
are c-number functionals which will be specified below. 
The quantities dN au (t) are known as Poisson increments. 
They are independent, random numbers which take on 
the possible values or I and satisfy the relation 

dN av (t)dN^(t) = 5 afi S vll dN a u(t). (16) 

Under the condition that dN au (t) = 1 for a particular a 
and v the other Poisson increments therefore vanish and, 
by virtue of the Eqs. i|14|) and (|15p. the state vectors then 
carry out the instantaneous jumps 

tpv — ► -iL av A a ^ v , Xv — ► M av B a Xv (17) 

The expectation values of the Poisson increments are 
given by 

E(dN av (t)) = T av dt. (18) 

This implies that dN av (t) — 1 with probability T av dt 
and, hence, the jumps (|17|l occur at a rate T au , which 



will also be determined below. If, on the other hand, all 
Poisson increments vanish we have d^v(f) — F v dt and 
dXv(t) — G v dt : which means that the state vectors follow 
the deterministic drift during dt. 

Our next step consists in deriving a stochastic equation 
for the random operator R(t) defined in Eq. JSJl, which 
will then lead to an equation of motion for the expec- 
tation value @. Employing the calculus of PDPs one 
finds 

dR = |d$i}<$ 3 | + |*i)(d$ 2 | + |d*i)<d$ 2 |. 

The third term on the right-hand side of this equation 
involves the products dN a idNp2 of the Poisson incre- 
ments, which vanish by virtue of Eq. i|ltj|) . This means 
that the state vectors and |$2(*)) evolve indepen- 

dently and that we may write 

dR= |d$i)<$ 2 | + |*i)(«f*2|. (19) 

With the help of the stochastic differential equations 1|12|) 
and l|13fl the state vector increments are found to be 

\d$ v ) = dip v ® Xu + ipv ® dxu + dijj v ® dx v 

= (F„dt + dJ u ) ®Xu + ^u® (G v dt + dK v ) 
+dJ v ® dK v . 

On using the structure of the jump terms (|14fl and i|15|) 
and relation l|16f) the third term may be written 

dJ v ® dK v = ^2 (-iL av A a - I) Vv 

a 

<g> (M av B a - I) X vdN au 
= -dJ v ® Xv 

+ ^ (-iL av A a - I) Vv 

a 

®M av B a x v dN av , 

which leads to 

|d$„) = F v dt® X v (20) 
+tp v ® {^G v dt - ^2 dN avX }j 
-%Y^ L a V M av (A a ^ v ) (g) (B a Xv)dN av . 

a 

This equation provides an exact relation for the stochas- 
tic increments \d<fr u ). To ensure that the first and the 
second term on the right-hand side vanish when taking 
the average over the Poisson increments, we now set 

F v = 0, G v = T vX u, (21) 

where 

r u = ^2T m , (22) 
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and 



1 



(23) 



This yields the expression 

|d*„) = i/) v ® f r v tft - ^ X. (24) 

a 

Finally, we substitute l|24ll into 1|19|) to arrive at 

= -i[Hi{t),R{t)]dt + dS(t). (25) 

Equation (|25Jl is the desired exact stochastic equation 
of motion of the random operator R(t) . The drift term in- 
volves the commutator with the interaction Hamiltonian 
Hj(t), while the noise term is given by the stochastic 
increment 



dS(t) = dTxR(t) + R(t)dTl 



(26) 



with 



dT v = P au dt - dN au ) (I + iT-lA a B a ) . (27) 

a 

According to Eqs. I|18|) and 127fl the average over the Pois- 
son increments yields E(<fT„) = 0. By virtue of Eq. I|26() 
this gives E(dS) = 0. Thus, if we take the average of both 
sides of Eq. (|2*5|l we are led directly to the von Neumann 
equation J3J. This shows that on average the stochastic 
dynamics defined by the differential equations (|12|) and 
(|13|l indeed reproduces the exact von Neumann dynam- 
ics of the density matrix of the combined system. We 
have thus achieved the goal of constructing a stochastic 
formulation of the evolution of the total system by means 
of a Markovian piecewise deterministic process. 

Up to this point the quantities L av and M av are com- 
pletely arbitrary with the only restriction that T au > 
(see Eq. which guarantees that the expectation val- 

ues ~Ei{dN av ) are positive, as it should be for random 
Poisson increments (see Eq. IjlHIl V In the following we 
choose 



WxA 



\B a xA 



(28) 



The advantage of this choice is that the jumps described 
by Eq. H±7|l then conserve the norm of the stochastic state 
vectors -0„ and Xv- Summarizing, the stochastic differ- 
ential equations defining the PDP now read as follows, 



\AM 



A a -I)ip v dN av (t), (29) 



dxu = T v Xvdt 



\B a Xu\ 



B a -I)x v dN av {t) 1 (30) 



where Y v is given by Eq. I|22|l and by 

WAcrtJW ■ \\ B *Xu\\ 



r 



^ll-IMI 



(31) 



We observe that ip v (t) is a pure, norm-conserving jump 
process, while Xv{t) is a PDP with norm-conserving 
jumps and a linear drift which leads to a monotonic in- 
crease of the norm of Xv 



B. Dynamics of fluctuations 

As a measure of the size of the fluctuations of the 
stochastic process constructed above we define 

D 2 {t) = E(\\R(t)-p(t)\\ 2 ) 

= E(tv{[R(t)-p(t)f[R(t)-p(t)]}).(32) 

The quantity D(t) is thus the root mean square dis- 
tance from the stochastic operator R(t) to its mean value 
pit) = E(i?(f)), the distance being determined through 
the Hilbert-Schmidt norm \ \A\ \ = ^/trjA^A}, where the 
trace is taken over the Hilbert space of the total system. 
Equation (|32|) may be written as 

D 2 (t) =E(tr{R\t)R{t)}) -trp 2 {t). (33) 

Since the dynamics of p(t) represents a unitary transfor- 
mation the trace over the square of p(t) is constant in 
time. For a pure initial state p(0) we have tip 2 = 1. 
Moreover, in the case of a sharp initial state, that is for 
R(Q) = p (0), one finds that D 2 {0) = 0. 

Our aim is to estimate the size of the fluctuations. 
To this end we first derive a differential equation for 
the mean square distance D 2 (t). With the help of the 
stochastic equation of motion (|25(l and of definition Ij26|l 
the differential of D 2 (t) is found to be 

dD 2 = E (tr {dS^dS}) 

= E(tr | dT}dTiRRt + dT%dT 2 R i R\y (34) 

Using then the definition l|27(l of the quantities dT v as 
well as Eqs. ©, dJ) and fl%)) . we obtain 



dD 2 „f^„ || (l + iT-lA a B a ) \<5> v )\\ 2 



e [J2r a J 



dt -\^-»» |||$ v )||2 

The choice <|2lfl) finally yields 



tr {R^R} 



dD 2 
dt 



2E ( ^r^trj^i?} 



(35) 



Equation l|35|) is an exact differential equation for the 
fluctuations of the random process. To find a rough esti- 
mate of the size of the fluctuations we suppose that the 
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rates T u are bounded from above, that is T u < Tq. This 
leads to the inequality 



dD 2 
~dT 



< 4r (D 2 + tvp 2 



which, on integrating, gives 



D 2 (t) < (trp 2 ) (e 



4r t 



L» 2 (0)e 



4r t 



(36) 



(37) 



This inequality provides a strict upper bound of the fluc- 
tuations of the random process. We note that the right- 
hand side of H37fl is finite for any finite time t. This leads 
to the important conclusion that the fluctuations of the 
process are finite for all finite times. 

Let us discuss in more detail the case of a sharp initial 
state, that is D 2 (0) = 0. We observe that for small times 
satisfying 4Tot <C 1 the root mean square distance then 
increases at most as the square root of time, 



D(t) < V(trp 2 )4IV. 



(38) 



For large times, AT t ^> 1, the root mean square distance 
may increase, however, exponentially with time, 



D(t) < Vtr^ 



(39) 



This shows that the stochastic method is useful for short 
and intermediate times, where the relevant time scale is 
given by l/2To. One further expects that the method 
is, in general, not efficient numerically for times which 
are large compared to l/2T , because of a possible ex- 
ponential increase of the fluctuations in this regime. It 
must be emphasized, however, that the statistical errors 
can be reduced considerably by employing the statistical 
independence of the increments |d$„) 



(see Sec. UTTSt . 
or by using a more complicated ansatz for the structure 
of the stochastic states (see Sec. 01. It should also be 
noted that the statistical errors are often much smaller 
than the upper bound given in the inequality <|39|) . An 
example will be discussed in Sec. IIVBI 



C. The stochastic simulation method 

1. Numerical algorithm 

The stochastic simulation method consists in a numer- 
ical Monte Carlo simulation of the stochastic differential 
equations and (|3l7|) . A realizations ipv(t), Xv(t) 01 
the process can be generated by means of the following 
algorithm. 

1. Suppose that the last jump into states Vv(t), Xv{t) 
occurred at some time t. In the case that t is the initial 
time t = 0, these states are taken to be the initial states 
which must be drawn from the probability distribution 
representing the initial density matrix through p(0) = 
E(R(0)). 

2. The next jump takes place at time t + r, where 
the r is a stochastic time step, the random waiting time, 



which is to be determined from the cumulative waiting 
time distribution function 



F(t) = 1 - exp (- dsT„ 



(40) 



A random number r following this distribution can be 
generated, for example, by drawing a uniform random 
number rj G (0, 1) and by solving the equation 



r\ = exp 



t+T 



dsT„(s) 



(41) 



for r. In between the previous and the next jump, that 
is within the time interval [t, t + T] the realization follows 
the deterministic drift which is given by 



Mt') = Mt), 



(42) 



Xu{t') = x*(t)exp\J t daT v (8)j, (43) 

where t < t' < t + t. 

3. Select a particular jump, that is select a particular 
value of the index a with probability 



T av {t + T) 

E a r Q ,(t + T)' 



(44) 



The corresponding jumps of the state vectors at time t+T 
then amount to the replacements 



IKVvll 
llx.ll 



^ctXv- 



(45) 
(46) 



Repeating these three steps until the desired final time 
tf is reached on obtains a realization Vv(t), X>(t) of the 
process over the whole time interval [0, tf]. An important 
feature of this algorithm is that it works with a random 
time step the size of which is adapted automatically by 
the algorithm: For large rates the time steps become 
small, while small rates lead to an enhancement of the 
time steps. For example, if T v is independent of time we 
simply have 



1 is 



(47) 



In the case of a time-dependent rate r„(t) it may well 
happen that the exponent in Eq. I|41|l is bounded from 
below and that, therefore, the exponential function con- 
verges to a finite value q > as r goes to infinity. For 
such a case one distinguishes two cases. For rj > q one 
determines t from Eq. (|41|) . while for rj < q one sets 
t = oo in which case there will be no further jumps. An 
example of this latter case will be shown in Sec. IIII Bl 

Finally we remark that for a numerical implementation 
of the simulation algorithm it might be more convenient 
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to employ a PDP with time-independent rates Y v . To 
this end one replaces the stochastic differential equations 
lUD and CIU by 



-iA n 



-I)^ v dN av (t), 



(48) 



/ b^d/V^f), (49) 



with an appropriate choice for constant rates r QJ ,. The 
advantage of this method is that the random waiting time 
is then always given by the simple expression (|47(l . The 
size of the statistical fluctuations, however, can depend 
considerably on the choice of the T alJ . 



2. Estimation of observables 

Suppose one has generated, by means of the algorithm 
described above, a sample consisting of J\f realizations of 
the process labeled by an index r, 

\K(t))=iP r u(t)®X r »(t), r=l,2,...,tf. (50) 
The quantum expectation value 

0(t) = tx{dp(t)} = E ((^ 2 (t)\d\^(t))) (51) 



of an observable O of the total system can then be esti- 
mated with the help of the ensemble average 



(52) 



In view of Eq. the reduced system's density matrix 
ps(t) is given through the ensemble mean 

^(*) = ^El^W)^(*)Kx5(*)lxI(*)>. (53) 



As emphasized already, the |$„(f)) evolve indepen- 
dently. Thus, if |$i(0)) and | < &2 (0)) are independent, 
as it is the case for a sharp initial value, for example, 
the processes |$i(f)) and |3>2(f)) are statistically inde- 
pendent. This implies that Eq. 1)51(1 can also be written 
in the following equivalent way, 



0(f) = (* 2 (t)|0|*i(*)>, 



(54) 



where ^(f)) = E(\$ v (t))). This suggests estimating the 
quantum expectation value (|51|l by means of the alterna- 
tive expression 



(55) 



Of course, the formulae l(52J) and 155|l lead to the same 
results in the limit of an infinite number of realizations. 



However, for a finite sample the statistical errors may 
differ considerably. 

To illustrate the difference between the statistical es- 
timates given by ((52(1 and H55|) , it suffices to consider the 
case O = \ip){Lp\, where \ip) may be any fixed state of 
the total system. We introduce the random quantities 
a = (i^l'f'i) and b = (cp\$> 2 ), as well as the corresponding 
realizations a r = and b r = (ipl^)- Equation 1(52(1 

can then be written as 



Oi 



Af 



a r . 



(56) 



The corresponding statistical error is provided by the ex- 
pression 



yVar(a)^ Var(a) + 2|E(a)| ^ 



where 



Var(o) =E(o*o)- |E(o)| s 



(58) 



is the variance of a, which is equal to the variance of b. 
On the other hand, Eq. 1(55(1 leads to the expression 



°2 = ^E^'- 



(59) 



The usage of this formula for the estimation of O is more 
efficient, in general, since the corresponding statistical 
error 



02 = 



/f?#WF (60) 



is smaller than u\ . The second method based on Eq. 1(59(1 
is thus to be preferred since it yields considerably smaller 
fluctuations. This difference between both methods be- 
comes particularly important if |E(a)| 2 , the quantity to 
be estimated, is small. The simulations presented in 
Sec. 1111 Bl and Sec. 1111 CI for example, have been carried 
out using this second method. 



3. Quantum correlation functions 

The fact that the stochastic method involves a pair 
of random wave functions also enables the design of an 
exact method for the determination of multitime corre- 
lation functions. The underlying idea is similar to the 
one employed in ^| for the calculation of correlation 
functions of quantum Markov processes. 

We restrict the discussion to the case of an arbitrary 
two-time correlation function of the form (X(t)Y(0)). In 
the interaction picture we can write (assuming t > 0) 



(X(t)Y(0)) = tr(X(t)U(t)Y(0)p(0)UHt)) 

= E(($ 2 (f)|x(f)i$r (*)», 



(61) 
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where X(t) and Y(t) are arbitrary operators in the in- 
teraction picture, and U(t) denotes the interaction pic- 
ture time-evolution operator of the total system over time 
t. The second line in Eq. (|61[) provides the stochas- 
tic representation of the quantum correlation function. 
In this expression both |$j (t)) and \$2(t)) follow the 
stochastic dynamics developed in Sec. Ill Al However, 
while the initial state of | *I»2 (*) ) is 1 3*2 (0)) , the stochas- 
tic process |<£>^(i)) evolves from the new initial state 
|$f(0)) = Y(0)|$i(0)). With this modification the 
stochastic algorithm for the determination of the corre- 
lation function is the same as above. The method can 
easily be generalized to the case of multitime correlation 
functions. An example will be studied in Sec. IIII Bl 



III. DECAY INTO A BOSONIC RESERVOIR 

To illustrate the general method developed in Sec. [n] 
we first study the model of a two-state system with ex- 
cited state |e), ground state \g), and corresponding tran- 
sition frequency ujq. This system is coupled to a bosonic 
reservoir consisting of field modes which will be labeled 
by an index k. The corresponding field operators that 
annihilate and create particles of frequency uik are de- 
noted by bk and b k , respectively. The interaction picture 
Hamiltonian is taken to be of the form 



Hi{t) = a + B(t)+<r-Bl(t) 



(62) 



The operators u+ = |e) (g\ and er_ = \g) (e| are the raising 
and lowering operators of the two-state system, while the 
reservoir operator B(t) is given by 



B{t) = y^fffc6 fc e 

k 



i(ui — uj k )t 



(63) 



with mode-dependent coupling constants gu- As a simple 
example we investigate the initial state 



|*„(0)>=V(0)®x(0) = |e>®|0>, 



(64) 



where |0) denotes the vacuum state of the reservoir. This 
initial state is statistically sharp and corresponds to the 
density matrix p(0) = |e)(e| <g> |0)(0| of the total system. 
This model can be solved analytically. The central phys- 
ical quantity that determines the influence of the reser- 
voir modes on the reduced system dynamics is provided 
by the bath correlation function 



f(t'-t) 



\B(t')B\t)\0) 
= / (LjJ{lo) exp[i(cjo — u){t' — £)], 



(65) 



which has been expressed here in terms of the spectral 
density J(ui). 



A. Description of the algorithm 

In the notation of Sec. Ill Al we have a 



1,2 and 

er_, Bi(t) = B{t) and B 2 (t) = B^(t). 



The application of the general technique of Sec. Ill C II 
to the present case leads to the following algorithm of 
simulating the stochastic dynamics. 

After an even number of jumps the reservoir state \v 
is proportional to the vacuum state. We thus infer from 
Eq. i|31|) that the transition rates are given by 



r*(f) 



\\B^t') X At)\\ 

\\Xu(t')\\ 



\B\t')\Q)\\ = v7(6)- (66) 



Since these rates are constant in time the random 
time step r is determined by Eq. (|47|l . that is r = 
— In 77/ yj /(0) with a uniform random number r\ in the in- 
terval (0, 1). Suppose that the previous jump took place 
at time t. Over the time interval [t, t + r] the state Xv 
then changes continuously according to 



Xv{t') = Xv {t)t 



t<t'<t + T, 



(67) 



until at time t - 
(|4^|l occur, 



t the jumps described in Eqs. I|45(l and 



-io--ip u (t + t), 
St(t + r) 



-xAt 



(68) 
(69) 



Note, in particular, that Xv jumps into a 1-particle state. 

After an odd number of jumps the reservoir state Xv 
represents a 1-particle state which was created out of the 
field vacuum at the time t of the last jump. Invoking 
again Eq. H31[) we find that the transition rates are now 
given by 



r„(f) = 



W) X v{t')\\ \\B(t')BHt)\0)\\ 



WxAt')\ 
\f(t'-t)\ 



v7(o) 



(70) 



We observe that these rates are time-dependent such that 
the random time step r as well as the deterministic drift 
of Xv must be determined from Eq. (|41|l and 1431) , respec- 
tively. In the present case we thus have 



X] = exp 



ds\f{a)\/>/W) 



(71) 



and 



rt'-t 



XV (*0 = X^)exp 



ds\m\/V7(fi) ■ (72) 



Finally, the jumps at time t + r take the form: 

ipu(t + r) — > -ia + tp u (t + T), (73) 

Xu(t + r) - 



B{t + T )xAt + T). (74) 



At time t+T the environment thus jumps back into a state 
which is proportional to the vacuum state. In terms of 
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Xv(t), which is defined to be the reservoir state just before 
the previous jump at time t, we can write the transition 
JHJ) as 



Xu(t + T) 



fir) 
\f{r)[ 



*„(t)exp / ds\f(s)\/Vf(0) 



This algorithm will be applied in the following two sec- 
tions to the damped Jaynes-Cummings model on reso- 
nance and with a finite detuning. 



B. Damped Jaynes-Cummings model on resonance 

The spectral density of the damped Jaynes-Cummings 
model on resonance is given by 



1 



7oA 2 



2tt (uj - oj) 2 + A 2 ' 



which yields the bath correlation function 



fit' - t) 



1 \ -\\t'-t\ 
-7oAe 1 i. 



(76) 



(77) 



This model can be used to describe the coupling of a two- 
level atom to an electromagnetic cavity mode which in 
turn is coupled to the continuum of modes of the electro- 
magnetic field vacuum. The quantity A -1 is the correla- 
tion time of the reservoir, while 7q can be interpreted 
as the Markovian relaxation time of the open system. 




FIG. 1: Excited state probability p(t) (Eq. of the 

damped Jaynes-Cummings model. Symbols: Monte Carlo 
simulations of the stochastic differential equations 12911 and 
with M = 5 ■ 10 6 realizations for the parameters A -1 = 
STo" 1 (diamonds) and A -1 = 207^ 1 (squares). The corre- 
sponding analytical solutions are given by the continuous and 
the broken line. 

The application of the simulation algorithm detailed 
in Sec. If 11 Al to this situation is straightforward. In par- 
ticular, we note that according to Eqs. (|70|l and l|77|l the 



waiting time distribution l|40(l after an odd number of 
jumps takes the form 



F(r) = l-expl-^[l-e-^]j. (78) 
Hence, the probability that no further jumps occur equals 

(79) 

This means that in the case 77 < q no further jumps 
occur, while in the case r\ > q the random time step is 
determined by Eq. (f7T|l which yields 

(80) 





0.5 



o 




-0.5. 



FIG. 2: The correlation function c(t) (Eq. (g2J) of the 
damped Jaynes-Cummings model: Analytical solution (con- 
tinuous line) and Monte Carlo simulation of the stochastic dif- 
ferential equations 1291 and 13011 (diamonds) for A^ 1 = 57 ( J~ 1 
and TV = 10 7 realizations. 

Results of Monte Carlo simulations of the damped 
Jaynes-Cummings model are presented in Fig. ^ which 
shows the population of the excited state, 



p(t)=E«e|^i)(^ 2 |e)<X2|xi» : 



(81) 



estimated from a sample of realizations of the stochastic 
process using the estimator described by Eq. (|59|l . As can 
be seen from the figure, the simulation results reproduce 
the analytical curves with high accuracy. We note that 
for the parameter values chosen the reservoir correlation 
time A -1 is larger than the reduced system's Markovian 
relaxation time 7o~ • We therefore observe a pronounced 
non-Markovian behavior and large deviations form the 
Born-Markov dynamics. For small and intermediate cou- 
plings, the open system dynamics derived from the model 
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described by the interaction Hamiltonian l|62[l and initial 
conditions (|64ll satisfies a time-local master equation of 
the form ps{t) — JC{t)ps(t) with a time-dependent super- 
operator JC{t). However, the TCL expansion of the gen- 
erator IC(t) breaks down in the strong coupling regime 
given by A -1 > ^l^ 1 for times t > to, where to denotes 
the first positive zero of p(t). Beyond the singularity 
at t = to the TCL expansion of the master equation 
is therefore not capable of describing the reduced sys- 
tem dynamics which develops a long memory time of the 
order to. However, as is exemplified in the figure, the 
stochastic simulation is seen to describe correctly the full 
non-Markovian behavior of the reduced system even in 
the strong coupling regime. 

To give an example of the simulation of correlation 
functions we investigate the quantity ((T_|_(t)<7_ (0)) which 
can be determined with the help of the method described 
in Sec. Ill C~3l Figure [21 shows the simulation results for 
the quantity 

c(t) = e-**V+(<)M0)>, (82) 
which again nicely fit the analytical curve. 



CL 




0.5 



5 y t 1 15 

'0 

FIG. 3: The excited state probability p(t) (Eq. JHB) of the 
damped Jaynes-Cummings model with detuning: Analytical 
solution (continuous line) and Monte Carlo simulation of the 
stochastic differential equations l|29|l and l|30|l (dots and er- 
rorbars) for A -1 = 57q~ , A = 70 and J\T = 10 7 realizations. 



C. Jaynes-Cummings model with detuning 

If the cavity mode is detuned from the atomic tran- 
sition frequency by an amount A the spectral density 
becomes 



1 



7oA 2 



2tt (w - A - lu) 2 + A 2 



which leads to the reservoir correlation function 
f(t'-t) = ^oAe^*'-*^!*'-*!. 



(83) 



(84) 



We can again use the simulation algorithm described in 
Sec. IIII Al although, by contrast to the previous case, the 
correlation function (|84|l is complex-valued. Since the 
transition rates and the deterministic drift of the process 
depend on the absolute value of /, the only modification 
of the algorithm for the resonant case appears in Eq. I|75|) 
which describes the even jumps into the vacuum state. 

An example of the simulation results is shown in Fig.|3 
The detuning A influences both the coherent dynamics 
of the system as well as the dissipation mechanism. This 
leads to a slower decay and to an oscillatory behavior 
of the excited state probability, which is correctly repro- 
duced by the stochastic simulation. 



IV. INTERACTION WITH A SPIN BATH 



specific central spin model which may be used to model 
the interaction of a single electron spin confined to a 
quantum dot with a bath of nuclear spins . 



A. Description of the model 

The model is defined by the total Hamiltonian 

N 



(85) 



The central spin is represented by the Pauli spin operator 
<t, while the N bath spins are given by the spin operators 
with j = 1,2,..., N. The coupling of the central 
spin to the jth bath spin is described by the constant 
. For simplicity, the coupling constants are taken to 
be — A/y/N. The corresponding interaction picture 
Hamiltonian can be written as 



Hi(t) = a 3 B 3 (t) + <J+B_{t) + <J_B + (t) (86) 



with 



B 3 = $>'V« 

j 

B ± = ^22A^a2 ] e^ ot . 



(87) 



The stochastic method developed in Sec. |H]is not re- 
stricted to the treatment of bosonic reservoirs. It is also 
applicable to the dynamics of open systems coupled to 
spin environments. As an example, we examine here a 



Our aim is to determine the coherence of the central 
spin, 



P+-(t) = (+\ps(t)\-), 



(89) 
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where |±) are the eigenstates of the 3-component 03 
of the central spin a with eigenvalues ±1. Within the 
stochastic simulation technique this quantity is repre- 
sented through the expectation value (see Eq. Q) 



p + _(i)=E«+|v> + )<iM-)<x-lx+>), 



(90) 



where we write here |<I>„) = |$±) = ip± ® x± f° r the 
stochastic states, that is the index v takes on the values 
v = i. The initial state is taken to be 



p(o) = l+X-l ® ^Ie- 



(91) 



denotes the unit matrix in the 2 Ar -dimensional state 
space He of the spin bath. The spin bath is thus in an 
unpolarized initial state. 



B. Simulation algorithm and results 

To apply the simulation technique it is useful to realize 
the unpolarized initial state 2~ N Ie of the spin bath with 
the help of an appropriate set of basis states of the Hilbert 
space He spanned by the N bath spins. To this end, we 
introduce states \ j, m) which are defined as simultaneous 
eigenstates of the square J 2 of the total spin angular 
momentum J of the bath and of its 3-component J3. 
The initial state can then be represented by 



|$±(0)> = \±)®\j,m) 



(92) 



with an appropriate probability distribution of the corre- 
sponding quantum numbers j and m which will be con- 
structed below. 

The state |$± (0)) defined in (|92|l is an eigenstate of the 
3-component ^(73 + J3 of the total spin angular momen- 
tum, which is a conserved quantity, corresponding to the 
eigenvalue |(±1 + 2m). This fact enables us to carry out 
the canonical transformation | $±(t)) — ► |$±(t)) defined 
by 



|*±(t)> =exp 



iAt 



((±1 + 2m)a 3 - 1) 



l*±(*)>, (93) 



which transforms the interaction Hamiltonian (|86[l into 

Hi(t) = a+B_ (t) + a- B+ (t) . (94) 

In this equation the B±(t) are given again by Eq. (|88ll . 
where, however, ujq must be replaced by the new frequen- 
cies u>±: 

2A 



UJQ 



(±l + 2m) 



(95) 



In terms of the stochastic states |$±) = rp± <8> x± the 
coherence of the central spin is then given by the expec- 
tation value 



P+-(t) = E (e^ Amt /^(+\^ + )^-\-)(X-\x + )) 



(96) 



Summarizing, we can simulate, employing the method 
developed in Sec.[nj the stochastic dynamics correspond- 
ing to the new interaction Hamiltonian (|94(l and estimate 
the coherence by means of the formula l|9t)l) . The canon- 
ical transformation l|93|) is accounted for in this formula 
by the exponential factor exp[— MAmt/y/N]. 

In order to see more explicitly how the method works 
it may be instructive at this point to consider first the 
simpler model obtained by omitting the terms cr±B^(t) 
of the interaction Hamiltonian (|86|l . The transformed 
Hamiltonian 194|) is then identically zero and the expres- 
sion H96JI for the coherence of the central spin becomes 



+N/2 _ 

E -AiAmt/VN 
Pm c 1 

m=-N/2 



(97) 

where p m is the probability of finding a basis state with 
quantum number m in the unpolarized initial mixture. 
Since all basis states are equally likely in this initial mix- 
ture, p m is found to be 



Pn 



2 N [± 



N 



Here, 2 N is the total number of basis states of the bath 
of N spins (the dimension of He), while the binomial co- 
efficient counts the number of basis states corresponding 
to a given value of m. The summation in Eq. I|97|l can 
easily be carried out to give 



P+-(t) 



2At 



N 



(99) 



which is the exact expression for the coherence of the 
central spin. We note that this expression may be ap- 
proximated by 



-2i4rf 



(100) 



in the limit of a large number of bath spins, N — ► 00, 
showing an exponential decay of the coherence of the 
central spin. Thus we see that the stochastic simulation 
for this simplified model reduces to the generation of a 
binomially distributed random number m and to the es- 
timation of the expectation value (|97|l . 

We turn again to the discussion of the full model de- 
scribed by the Hamiltonian (|86|l . Employing the method 
described above and using the transformed interaction 
Hamiltonian l|94l) we see that the simulation algorithm 
is quite similar to the one used already in the bosonic 
case. In fact, the simulation technique turns out to be 
even simpler. Suppose we have drawn the initial state 
|±) (g) \j, m). The bath state X±{t) then jumps between 
states which are proportional to \ j, m) and \ j, mil). The 
corresponding jump rate 



2A 



j(j + 1) — m(m i 1) 



N 



(101) 
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is independent of time. The waiting time of the PDP is 
therefore always exponentially distributed, which makes 
the numerical implementation particularly easy for this 
case. A detailed analysis of the process reveals that the 
coherence can be represented through the expectation 
value 

p+-(t) = E(e~ 4lAmt/VW (-l) ( ' £++fe - )/2 e (r + +r - )t 
x exp(^ + (r+ +t+ ... + r+ + )) 
x exp(wj_ (t 2 - + t 4 - . . . + r,7_ ))) . (102) 

Here, denotes the random time step before the 2nth 
jump of |$±), while r± and u± have already been defined 
in Eqs. and Ij95|l . The quantity k± is defined as the 

total number of jumps of \<$>± (t)) during the time interval 
from to t. The integers k± may be supposed to be even 
since only trajectories with an even number of jumps 
contribute to the expectation value (|1U2|) . 

It remains to explain how to generate, in the general 
case, the initial states \j, m) in Eq. 1(52)1. More precisely, 
these states should be written as |A, j, m), where A stands 
for an additional quantum number which, together with 
j and to, uniquely fixes the basis state. The quantum 
number A corresponds to further observables of the spin 
bath which commute with J 2 and J3. If N is even j takes 
on the values j = 0, 1, 2, . . . , ^, while j = i, |, . . . , =j if 
N is odd. For a given value of j the quantum number m 
takes on the values to = — j, + 

In order to achieve that the initial ensemble represents 
the unpolarized bath state, that is 

E(\X,j,m){X,j,m\) = ^I E , (103) 

all basis states | A, j, to) must occur with the same proba- 
bility of 2~ N . Since the value of the quantum number A 
is irrelevant in the simulation scheme, we need the prob- 
ability P(J,m) of finding the pair of quantum numbers 
(j, m) in the initial ensemble. This probability can be 
written as 

P(j,m) = 2- N af . (104) 

The quantity denotes the number of times a given an- 
gular momentum j appears in the decomposition of the 
Hilbcrt space He of N spins into irreducible subspaces of 
the rotation group. Since a certain j-manifold consists of 
(2j + 1) states, distinguished by their values of the quan- 
tum number to, we can also say that (2j + l)a^ is equal 
to the number of independent ways the TV bath spins can 
be coupled to give the total angular momentum j. For 
example, the Hilbert space of N — 4 spins decomposes 
into two (j = 0)-manifolds, three (j = l)-manifolds, and 
one (j = 2)-manifold, that is we have ag = 2, a\ = 3, 
and af = 1. It may be shown |20| that is given by 
the general expression 



We note that P(J, m) is normalized, 
+3 

J2 E p (i' TO ) = 1 . ( 106 ) 

j m=—j 

and does of course not depend on to. In summary, the 
quantum numbers (j, m) of the initial ensemble follow 
the distribution P(j, to) given by the expressions l|104|) 
and (|105fl . In the stochastic simulation algorithm one 
therefore has to generate a sample of random numbers 
(j, to) with this distribution, which is easily done making 
use of the inversion method, for example. 
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FIG. 4: Real part of the coherence l|89|l of the central spin in- 
teracting with a spin bath through the Hamiltonian 186II with 
N = 10 3 . Symbols: Monte Carlo simulation of the stochas- 
tic differential equations 1291 and 130H using Af = 2 ■ 10 7 
realizations for the parameters A/luo = 0.1 (diamonds), 
A/ujo ~ 0.2 (squares), and A/u>o = 10 (triangles). Continuous 
lines: Corresponding solutions of the von Neumann equation 
©. The dashed line (A/uio = 0.1), the dashed-dotted line 
(A/ujo = 0.2), and the dotted line (A/u>o — 10) show the re- 
sults obtained from the TCL master equation in second order 
(Eqs. and (HEJ). 

Examples of Monte Carlo simulations of the central 
spin model are shown in Fig. ^ One observes that the 
PDP reproduces the von Neumann dynamics with high 
accuracy. We do not show errorbars in the figure be- 
cause the statistical errors are smaller than the size of 
the symbols. The figure also displays the results found 
with the help of the second-order TCL master equation 
of the central spin which is given by 

d „ , 9 1 — cos Wnt . . , „„. 

-ps = -2iA 2 -[<J 3 , P s] (107) 

at ujq 

~A 2 t[a 3 ,[cr 3 ,ps}} 

. . 2 sin wo i / l r t 

+4A a-psv+ - -{o+(J-,ps\ 

uj \ 2 

+a + p s a- - )-{(T-a +1 p s } J . 
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The solution of this master equation is easily constructed. 
It yields the expression 

P+ _(t) = cxp[-r(<)] P+ _(0) (108) 

for the coherence of the central spin, where 

= AJJHS _smu*A 

uj q \ UJ Q t J 

For the parameter values chosen the exact dynamics of 
the central spin is seen to deviate significantly from the 
one predicted by the second-order TCL master equation. 
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FIG. 5: Statistical errors a(t) of Monte Carlo simulations of 
the central spin model with 10 7 realizations, A/uiq = 0.5 and 
three different values of the number of bath spins: N = 10 
(diamonds), N = 100 (squares), and N = 1000 (triangles). 
The continuous line shows the estimate given by Eq. Ill 11 . 

Figure presents an example of the behavior of the 
fluctuations of the stochastic process. The figure shows 
a plot of the statistical errors a(t) of three Monte Carlo 
simulations with a fixed number Af of realizations, but 
with three different values of the number N of bath spins. 
We conclude from the figure that, within the range of 
time investigated, a(t) is roughly independent of N. To 
understand this behavior we refer to expression l|102|) 
which yields 

/E(cxpf2(r+ + r_)i) , . 

<r{t) < y ^ P[ y ^. (110) 

The right-hand side of this inequality may be estimated 
by replacing the random quantities r± by suitable aver- 
ages using the distribution l|104|) . This gives the estimate 

„( t )~22j^I. (in) 



This expression is indeed independent of N and provides 
a good estimate of the standard error in the given time 
interval, as can be seen from the figure. Moreover, this re- 
sult implies that the fluctuations grow with a rate which 
is much smaller than the one provided by the strict upper 
bound 2Fo of T + In fact, Tq scales with the square 

root of N which predicts a much stronger increase of the 
fluctuations. 



V. CONCLUSIONS 

It has been shown in this paper that the von Neumann 
dynamics of a combined quantum system can be formu- 
lated in terms of a rather simple piecewise deterministic 
process which gives rise to a powerful and efficient Monte 
Carlo simulation method of the exact non-Markovian re- 
duced system behavior. A Markovian representation of 
the dynamics was achieved through the use of a pair of 
product states i/j v ® \v m t he state space of the total sys- 
tem. The stochastic propagation of an ensemble of such 
pairs then enables one to mimic the exact time-evolution 
of the reduced system's density matrix. 

The examples discussed in Sec. HHI and llVl illustrate the 
generality of the method: It is applicable to both bosonic 
and spin environments and is not restricted to linear dis- 
sipation or to a perturbation treatment of the system- 
environment coupling. Most importantly, the method 
does not require the derivation, not even the existence of 
a master equation of the reduced system. At the same 
time, the technique allows the direct determination of all 
kinds of multitime quantum correlation functions. Al- 
though our discussion was carried out in the interaction 
picture, it is obvious that the stochastic dynamics can 
also be formulated in the Schrodinger picture, in which 
case both ?/v and Xv follow, in general, a non-trivial de- 
terministic evolution. Furthermore, it should be clear 
that, instead of using a PDP, one can also employ a dif- 
fusion process (Brownian motion) to construct an unrav- 
eling of the von Neumann equation. 

The stochastic technique was formulated here as a 
method of simulating the dynamics of open systems in 
real time. A potential extension of the method is to re- 
formulate the dynamics in imaginary time |2lj . in or- 
der to determine the properties of the system in ther- 
modynamic equilibrium. With the total Hamiltonian 
H = Hs + He + Hi in the Schrodinger picture the canon- 
ical equilibrium density matrix (not normalized) is given 
by p(/3) — e~P H , where (i — l/k^T is the inverse tem- 
perature. At infinite temperature we have p((3 = 0) = I. 
This suggests determining the equilibrium density at fi- 
nite temperature by solving the evolution equation 

± p ( s ) = -^{H,p(s)} (112) 

over the interval from s = to s — (3. This imaginary- 
time dynamics can again be represented in terms of a 
stochastic process for a pair of product states \&„(s)) = 
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® Xv( s )- An appropriate system of stochastic dif- 
ferential equations in the Schrodinger picture is given by 

+ (-\ L ^A a - J J i> v dN av , (113) 
+ ^ {M av B a - /) XudN av . (114) 

a 

Performing a calculation analogous to the one of Sec. Ill Al 
it is easy to verify that the expectation value p(s) = 
E(|$i(s))($2(s)|) satisfies the evolution equation (|112[) . 
The dN a „(s) are again independent Poisson increments 
satisfying ~E(dN av (s)) — T al/ ds, and the relations i|22|) 
and 

remain valid. 

An important restriction of the Monte Carlo technique 
is provided by the behavior of the statistical fluctuations. 
The considerations of Sec. Ill Bi as well as the example dis- 
cussed in Sec. lIVBl reveal that the method as formulated 
in Sec. Ill Al is feasible, in general, only for short and in- 
termediate time scales. For large times statistical errors 
may grow exponentially fast, ruling out the estimation 
of statistical quantities with reasonable effort. However, 
this conclusion rests on the assumption that the stochas- 
tic states |$„(i)) are tensor products of certain system 
and environment states. This leads to a further poten- 
tial generalization of the method, namely to introduce a 
class of stochastic states with a more complicated struc- 
ture, the aim being a more efficient representation of p(t) 
as the expectation value over the corresponding random 
process. 

Since the interaction generally creates correlations be- 
tween the states of system and environment it could be 
advantageous, e. g., to use a class of entangled stochastic 
states. The spin bath model studied in Sec. |^ leads to 
a trivial example: The class of entangled states defined 
by (a and (3 are complex amplitudes) 

a\+)®\j,m)+l3\-)®\j,m + l) (115) 



yields an extremely efficient stochastic representation of 
the dynamics: As a consequence of the conservation of 
the 3-component of the total spin angular momentum, 
the subspaces spanned by the states |+) ® \j, m) and 
|— ?7i+l) are invariant under the time-evolution and, 
thus, the dynamics may be expressed entirely though an 
appropriate (deterministic) time-dependence of the am- 
plitudes a and (3. Therefore, only the initial state is a 
random quantity and the statistical errors are constant 
in time. 

In a further possible extension of the method one could 
employ a stochastic evolution of mixed states instead of 
pure states. As an example we introduce a stochastic 
matrix 

fl(t) = HM*)}<$i(*)l® (us) 

where the ipi>(t) are random states of the open system 
and Rs(t) is a random operator in He, and try again to 
find stochastic evolution equations such that the exact 
von Neumann dynamics is recovered by means of the ex- 
pectation value p{t) = E(i?(t)). This is indeed possible 
if we use the stochastic differential equations 129fl for the 
i/j v {t) and if we replace l|3(J|) by the following stochastic 
differential equation for the random operator Re it), 

dR E = TR E dt + Y {M al B a - I) R E dN al 

a 

+ J2 R e (M a2 Bi - I) dN a2 , (117) 

a 

where T = Ti+r2 = Ylav ^ av 1S ^ ne total jump rate. The 
further development of the stochastic technique proposed 
in this paper should include a systematic investigation of 
the potentialities of the extensions indicated above. 



Acknowledgments 

The author would like to thank F. Petruccione for help- 
ful discussions and comments. 



[1] C. Cohen- Tannoudji, J. Dupont-Roc, and G. Gryn- 
berg, Atom-Photon Interactions (John Wiley, New York, 
1998). 

[2] C. W. Gardiner and P. Zoller, Quantum Noise (Springer- 

Verlag, Berlin, 2000). 
[3] S. Nakajima, Progr. Theor. Phys. 20, 948 (1958). 
[4] R. Zwanzig, J. Chem. Phys. 33, 1338 (1960). 
[5] R. P. Feynman and F. L. Vernon, Ann. Phys. (N.Y.) 24, 

118 (1963). 

[6] A. O. Caldeira and A. J. Leggett, Physica 121A, 587 
(1983). 

[7] H. P. Breuer and F. Petruccione, The Theory of Open 



Quantum Systems (Oxford University Press, Oxford, 
2002). 

[8] G. M. Moy, J. J. Hope, and C. M. Savage, Phys. Rev. A 
59, 667 (1999); Phys. Rev. A 61, 023603 (2000). 

[9] G. M. Palma, K.-A. Suominen, A. K. Ekert, Proc. R. 
Soc. Lond. A 452, 567 (1996). 
[10] N. V. Prokof'ev and P. C. E. Stamp, Rep. Prog. Phys. 
63, 669 (2000). 

[11] J. Dalibard, Y. Castin, and K. M0lmer, Phys. Rev. 
Lett. 68, 580 (1992); R. Dum, P. Zoller, and H. Ritsch, 
Phys. Rev. A 45, 4879 (1992); N. Gisin and I. C. Per- 
cival, J. Math. Phys. A: Math. Gen. 25, 5677 (1992); 



14 



H. Carmichael, An Open Systems Approach to Quantum 

Optics, Lecture Notes in Physics ml8 (Springer- Verlag, 

Berlin, 1993); M. B. Plenio and P. L. Knight, Rev. Mod. 

Phys. 70, 101 (1998); A. Imamoglu and Y. Yamamoto, 

Phys. Lett. A 191, 425 (1994); A. Imamoglu, Phys. Rev. 

A 50, 3650 (1994). 
[12] L. Diosi, N. Gisin, and W. T. Strunz, Phys. Rev. A 58, 

1699 (1998); W. T. Strunz, L. Diosi, and N. Gisin, Phys. 

Rev. Lett. 82, 1801 (1999). 
[13] H. P. Breuer, B. Kappler, and F. Petruccione, Phys. Rev. 

A 59, 1633 (1999). 
[14] J. T. Stockburger and H. Grabert, Phys. Rev. Lett. 88, 

170407 (2002). 



[15] H. P. Breuer, quant-ph/0308052. 

[16] I. Carusotto, Y. Castin, and J. Dalibard, Phys. Rev. A 

63, 023606 (2001). 
[17] O. Juillet and Ph. Chomaz, Phys. Rev. Lett. 88, 142503 

(2002). 

[18] H. P. Breuer, B. Kappler, and F. Petruccione, Phys. Rev. 

A 56, 2334 (1997). 
[19] A. V. Khaetskii, D. Loss, and L. Glazman, Phys. Rev. 

Lett. 88, 186802 (2002). 
[20] A. Hutton and S. Bose, quant-ph/0208114. 
[21] I. Carusotto and Y. Castin, J. Phys. B: At. Mol. Opt. 

Phys. 34, 4589 (2001). 



